An extended analysis of the viscosity kernel for monatomic and diatomic fluids 
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We present an extended analysis of the wave-vector dependent shear viscosity of monatomic and 
diatomic (liquid chlorine) fluids over a wide range of wave-vectors and for a variety of state points. 
The analysis is based on equilibrium molecular dynamics simulations, which involves the evaluation 
of transverse momentum density and shear stress autocorrelation functions. For liquid chlorine we 
present the results in both atomic and molecular formalisms. We find that the viscosity kernel 
of chlorine is statistically indistinguishable with respect to atomic and molecular formalisms. The 
results further suggest that the real space viscosity kernels of monatomic and diatomic fluids depends 
sensitively on the density, the potential energy function and the choice of fitting function in reciprocal 
space. It is also shown that the reciprocal space shear viscosity data can be fitted to two different 
simple functional forms over the entire density, temperature and wave-vector range: a function 
composed of n-Gaussian terms and a Lorentzian type function. Overall, the real space viscosity 
kernel has a width of 3 to 6 atomic diameters which means that the generalized hydrodynamic 
constitutive relation is required for fluids with strain rates that vary nonlinearly over distances of 
the order of atomic dimensions. 



I. INTRODUCTION 



a local viscosity defined by Newton's viscosity law as 
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Fluid dynamics at atomic and molecular scales still 
present a challenge for theoreticians as well as for experi- 
mentalists. Molecular dynamics (MD) is a computational 
tool that has contributed significantly to the fundamen- 
tal understanding of these systems by providing informa- 
tion about processes not directly approachable by exper- 
imental studies. A central problem in the study of fluids 
at such small length and time scales is the computation 
of meaningful transport properties. Many equilibrium 
molecular dynamics (EMD) as well as nonequilibrium 
molecular dynamics (NEMD) simulations of nanofhiids 
have been performed since the early 1980s [H-||. In most 
of these simulations the stress was treated as being de- 
pendent of the local strain rate rather than the entire 
strain rate distribution in the system. Todd et al. have 
recently shown that in all but the simplest flows (e.g. pla- 
nar Couette and Poiseuille flows) and for velocity fields 
with high gradients in the strain rate over the width of 
the real space viscosity kernel, non-locality can play a 
significant role [1, 0] . In the case of a homogeneous fluid, 



P xy (r,t) = - Vo 5(r-r',t-t')j(v',t')dr'dt' (1) 
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even exhibits singularities at points where the strain rate 
is zero (sT-flll]. I n Eq. ([T]) P xy (r,t) represents the (x,y) 
off-diagonal component of the pressure tensor, 7(1", t) is 
the shear strain rate at position r and time t, and rjo is the 
local shear viscosity. In general, a nonlocal constitutive 
equation that allows for spatial and temporal non-locality 
can be expressed as [H, HH 

t 00 

P xy (r,t) = -f J n(r-r',t-t')i(r',t')dr'dt', (2) 

-00 

for a homogeneous fluid. In the situation where the strain 
rate is constant in time and only varies with respect to 
the spatial coordinate y, Eq. ^ can be written as 



p xy(y) 



v(v - y')i(y')dy' 



(3) 
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In reciprocal space Eq. (j3]), can be expressed as 

Pxy(ky) = -fj(ky)^{ky), (4) 
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where k y is the y component of the wave-vector as de- 
fined later in Section IIIII Such a constitutive equation is 
expected to be necessary for the description of flows in 
highly confined systems, due to the large change in the 
strain rate with position in the vicinity of the wall 0. 

The best available theoretical predictions of the wave- 
vector dependent viscosity are based on mode-coupling 
theory and generalized Enskog theory (l^4l6l|. However, 
these theories do not quantitatively agree with data ob- 
tained via computer simulations [12fl . The theoretical 
predictions focus on the transverse momentum density 
autocorrelation function, which is found by an iterative 
numerical solution of a system of nonlinear equations. 
Consequently, the theories do not result in analytical ex- 
pressions for the correlation functions or the wave- vector 
dependent transport coefficients, which are the focus of 
the present study. More recently, a modified collective 
mode approach has been successfully applied by Omelyan 
et al. [ljj to the TIP4 model of water. In contrast 
to other semi-phcnomcnological approaches used for in- 
stance in TIP4P and SPC/E models of water by Bertolini 
et al. [U[ and Palmer [l9j|, Omelyan et al. reproduced 
the reciprocal space kernel using a relatively small num- 
ber of modes. 

In this paper, we extend the work done by Hansen et 
al. [IFJ] and focus on computing the spatially non-local 
viscosity kernel for monatomic and diatomic fluids over 
a wider range of wave- vectors, state points and potential 
energy functions. We arc specifically interested in identi- 
fying functional forms that fit the reciprocal space kernel 
data. On the basis of these results, we are be able to 
assess the length scale (i.e. the width of the real space 
kernel) over which the governing generalized constitutive 
relation Eq. must be used. We expect that non-local 
transport phenomena would be relevant in shock waves 
12 , l2ll - l2~i ] , shear banding [25| , flows of micellar solutions 
26], suspensions of rigid fibers [13] and jammed or glassy 
systems [28| . 

This paper is structured as follows: In Section III Al 
we give an overview of the general formulation and the 
expressions for the complex wave-vector and frequency 
dependent viscosity are given. In Section fllll we describe 
the simulation methodology and conditions. In Section 
IVI we present our molecular dynamics simulation data 
and compare the results of our monatomic and diatomic 
viscosity kernels, particularly the shape of the kernels. 
Finally, we summarize and conclude our analysis in Sec- 
tion 5. 



II. METHODOLOGY 

We shall briefly introduce the main conceptual back- 
ground used in this work, namely, the wave-vector de- 
pendent momentum density, stress and viscosity in the 
atomic and molecular formulations. 



A. Wave-vector dependent momentum density for 
atomic and molecular fluids 

For a single component atomic fluid the real space mi- 
croscopic momentum density is given by fl3j 



J(r,i) = p Q (r,i)v(r,i) 



N a 

E 



77ijVj(i)(5(r - Vi) (5) 



where p a (r,t) — m i^( r — r ») ^ s the mass density, 

m,, r, and Vj are the mass, position and velocity of atom 
i. The summation runs over the number of atoms N a 
in the system. The Fourier transform of the momentum 
density is 



J(k,t) = y^mjVj(f)e 

i=l 



(6) 



while the Fourier transform of the mass density is 



Pa(k,t) 



Ei=l m i e 



We define the Fourier 



transform of a function f(r) as T[f(r)] = /(fc) 



I. 
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f(r)dr. 



The atomic representation of the momentum de nsity 
for a molecular fluid can be written in real space as [29 ]: 



3 A {r,t) = p a (r,t)v(r,t) = X] E m iaV ia (t)5(r - r la ) 

i—l a — 1 

(7) 

where the mass density is defined as p a (r,t) = 

Si=a m iu&{ Y ~ r ia)- The inner summation extends 

over the N s mass points in a molecule and the outer sum- 
mation extends over the number of molecules N m in the 
system. In general, N s depends on the molecule index i 
for a multicomponcnt system, but in our systems N s is 
the same for all molecules and the interaction sites are 
assumed to have identical mass, namely m !Q . 

The Fourier transform of the momentum density is 



N m N 3 



(8) 



where the transformed mass density is p a (k,t) 



Etd Eq=i ™iae !k ' r,c * ■ For molecules composed of N s 
atoms we can define the mass of molecule i, Mi — 
53 *i m iai position of the molecular center of mass as 
rj = Yla=i m ia r ia/Mi, position of site a of molecule i 
relative to the center of mass of molecule i as Ri a = 
Ti a — Vi, and center of mass momentum of the molecule 

This means that the atomic 



as pi 
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mass density can be written in k-space as /5 a (k, t) = 

Et™ Eq=i «ii a e ik '' r * +R '°'- If we expand this rela- 
tion further we can express the atomic mass density 
in terms of molecular mass density in which we de- 
fine the mass density in the molecular representation 
as p m (k,t) = 53^J Mie lk ri in reciprocal space and 
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as p m (r,t) — ^Zi^i ^i<5(r — in real space, respec- 
tively. In a similar way we can expand the atomic mo- 
mentum density about the molecular center of mass: 

j A (k, t) = £«=i m ia v ia (i + zk • R iQ + . . . ) e « k -*° . 

The Fourier transform of the momentum density in the 
molecular representation can then be defined as 



N„ 



J M (M) = 5>/ lVl (i) ( 



(9) 



i=l 



A complete procedure for expressing the mass and mo- 
mentum densities in physical and reciprocal space for 
atomic and molecular fluids has been discussed in more 
detail by Todd and Daivis 



B. Wave- vector dependent pressure tensor 

For a monatomic system the wave-vector dependent 
pressure tensor is defined as [29| 

TV N N 

2—1 Z— 1 

(10) 

where Fj is the force on atom i due to atom j and <?(k) = 
(e lk r - l)/ik • r = J2n=o( ik ■ r) n /{n + 1)! is the Fourier 
transform of the Irving-Kirkwood Oij operator [30l ] . 

For a molecular system the molecular pressure tensor 
is the pressure calculated using the intermolecular forces 
and the molecular center of mass momenta. The atomic 
pressure on the other hand includes all atomic momenta 
and all interatomic forces and constraint forces. Thus the 
wave- vector dependent pressure tensor for constrained di- 
atomic fluid can be written in atomic representation as 

N m 2 



%—\ a—1 

^ N m 2 Nm 2 



i—1 a — 1 /3 — 1 
N m 2 

i=l a=l 



(11) 



where F ia jp is the LJ force acting on site a of molecule i 
due to site /3 of molecule j, and Ff a is the bond constraint 
force on site a of molecule i. = Tjp — Yi a is the 

minimum image separation of site a of molecule i from 
site (3 of molecule j. 

The pressure tensor for a diatomic molecule in the 
molecular representation is defined as 



N, 



P M (M) = £ 



Pi Pi ik-r,- _ 1 

Mi 6 2 

i—1 i—1 jt^v 



^^ rij r^(k)e fa fl2) 



where, F^ represents the total intermolecular force on 
molecule i due to molecule j. = Vj — r. ; is the mini- 
mum image separation of the center of mass of molecule 



i from the center of mass of molecule j. In cases where 
two sites on two different periodic images of the same 
molecule may interact, the correct value of = Tj — 
corresponding to the particular images of molecule i and 
j in Fiaj/3 must be used. Though this is unlikely to hap- 
pen in diatomic fluids it is particularly important in sim- 
ulation of long molecules. The momenta appearing in 
these equations, Pi Q , Pi, are the momenta appearing in 
the respective atomic and molecular equations of motion 
Eqs. (ETJ and CT. 



C. Wave- vector and frequency dependent viscosity 

The complex wave- vector and frequency dependent vis- 
cosity can be evaluated by using two different expres- 
sions in terms of the Fourier-Laplace transform of the 
transverse momentum density autocorrelation function 
(ACF), Cj_(k, t), and the Fourier-Laplace transform of 
the stress tensor autocorrelation function, 7V(k, t) fl3l ]. 
We define the complex Laplace transform (one-sided 
Fourier transform) as C[f(t)] = f(w) = J °° f(t)e~ lut dt. 
We also note that for the sake of simplicity of notation 
and consistency with the notation used in previous pub- 
lications, in what follows, we drop the tilde sign over the 
Fourier transformed correlation functions and keep the 
tilde notation over the Fourier-Laplace transformed cor- 
relation functions only. If we set = (0, k y ,0) and J x 
is the component of the momentum density in the x di- 
rection, the expression for the wave- vector and frequency 
dependent viscosity in terms of C±(k y ,uj) takes the form 

m 

, \ _ P C±(ky,t = 0) - iU]C±(ky,u) 



C±(k y ,u}) 



where p is the number density of the fluid and Cj_(k yi to) 
is the Laplace transform of the ensemble averaged 
transverse momentum density autocorrelation function 
C±(k y ,t), which is defined as 

C X (ky,t)= ~(j X (ky,t)J X (ky,t = 0)y (14) 

The zero time value of C±(k y ,t = 0) in the thermody- 
namic limit is 



C x (k y ,t = 0) =pk B T. 



(15) 



fcg is Boltzmann's constant. Due to finite time aver- 
aging and finite system size, the value of C(k y ,t = 0) 
obtained from simulation differs insubstantially from the 
exact value given by 



C±(ky,t = 0) =pk B T 



3N -4 
3N 



(16) 



because the total peculiar kinetic energy and three com- 
ponents of the momenta are constants of the motion in 
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our simulations. We also note that the number of degrees 
of freedom in the simulated system has not been taken 
into account in Eq. (fl~5|) . therefore we use the simulated 
value in our calculations to ensure numerical consistency 
of the computed properties. 

The expression for the wave-vector and frequency de- 
pendent viscosity in terms of the autocorrelation function 
of the shear stress N(k y ,u>) takes the form: 



fj{k y ,w) 

where 

N(k y ,uj) 



N{ky,Uj) 



C(k y ,t = 0)/pk B T - k 2 N(k y ,uj)/iujp 



(17) 



Vk B T 



X 



P Xy (ky,t)P X y(k y ,0)) . (l8) 



In the zero wave- vector limit, a generalization of the 
Green-Kubo expression for the shear viscosity allows the 
transverse momentum flux to be in an arbitrary direction 
rather than along a coordinate axis and can be written 
in terms of the stress tensor as [3l], H3] : 



V 



V 



10k B T 



P os {t) : P os (0) )dt 



(19) 



where the os superscript denotes the traceless symmet- 
ric part of the stress tensor P os (t) = ±[P(i) + P T (t)} - 
|fr[P(i)]l and V is the simulation volume. In an 
isotropic fluid, because the tensor, P os , appearing in 
Eq. (fT9"|) is traceless and symmetric, the shear viscosity 
is also traceless and symmetric. Consequently, the shear 
viscosity can be expressed in terms of the invariant / of 
the shear viscosity tensor as rj = 7/10. 

In the case u — > and k y — > 0, relation (fT?) reduces 
to the Green-Kubo formula [3l|. All the non-zero wave- 
vector integrals, Eq. (fT8|) . converge to zero [l3[ there- 
fore the relation in Eq. (|13|) must be used when non- 
zero wave-vector viscosities are calculated. By comput- 
ing the integrals one can computationally verify the zero 
values of the zero-frequency limit of the function N(k, lu) 
and thus demonstrate why neither substitution of uj = 
into Eq. (fT?} nor evaluation of Eq. (|T5)) at non-zero wave- 
vector yields the zero-frequency wave-vector dependent 
viscosity. 



III. SIMULATION DETAILS 

We use the Edberg, Evans, and Morriss algorithm 
[33M35I] with an improved cell neighbour list construc- 
tion algorithm [36] to perform equilibrium simulations 
at constant N,V,Tm- N is cither the number of atoms 
or molecules and Tm is the molecular temperature as 
defined by Eq. (f2"S"j) . For the atomic fluids studied in 
this work, the atoms interact via a Lennard- Jones (LJ) 
or Weeks-Chandler-Andersen (WCA) potential energy 
function (37| . The LJ atoms have an interaction poten- 
tial truncated at r c = 2.5a and WCA atoms have an 



interaction potential truncated at r c 
general: 
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2 1 / f V 



In 



.0- 



ij ^ Tc 
ij ^ 1 c 

(20) 

where is the interatomic separation, e is the poten- 
tial well depth, and a is the value of r%j at which the 
unshifted potential is zero. The shift $ c is the value of 
the unshifted potential at the cutoff r.y = r c , and is in- 
troduced to eliminate the discontinuity in the potential 
energy. At distances greater than the cutoff distance r c , 
the potential is zero. 

Our diatomic model of liquid chlorine is similar the 
the one used by Edberg et al. [38[, Hounkonnou et al. 
[39| , Travis et al. [10, [4l[ and more recently by Matin 
et al. [ID, HH to allow a direct comparison of our results 
with previous work. This model represents chlorine as a 
diatomic LJ molecule with r c = 2.5a and a fixed bond 
length of 0.63cr. For an adequate representation of the 
properties of chlorine the LJ parameters are: a = 3.332A 
and e/k B = 178. 3K. Liquid chlorine systems of 108 and 
864 molecules are studied at a reduced site number den- 
sity of p a = 1.088 and a reduced molecular temperature 
Tm = 0.970. We summarize the most important simu- 
lation parameters in Table [U All our simulations were 

TABLE I: Simulation details 





WCA 


LJ 


Chlorine 


Site number density,p n 


0.375,0.480,0.840 


0.840 


1.088 


Temperature, T 


0.765, f.O 


0.765, 1.0 


0.97 


Number of atoms, N a 


108,2048,6912 


108, 2048, 6912 


216, 1728 


Number of sites, N s 


1 


1 


2 


Bond length, / 






0.63 


LJ cutoff, r c 


2 i/c 


2.5 


2.5 



carried out in a cubic box with periodic boundary condi- 
tions. The fifth-order Gear predictor corrector algorithm 
[ii], |45[ with time step 5t = 0.001 was employed to solve 
the equations of motion. The equations of motion can be 
written for a monatomic fluid in the isokinetic ensemble 
(at equilibrium) as [l3[: 



El 



Pi - (APi 



(21) 



where i denotes atom i. is the position, is the 
momentum and m, is the mass of the designated atom. 
Fj is the force on atom i due to other atoms and £4 is the 
atomic thermostat multiplier. The thermostat multiplier 
is chosen so as to fix the kinetic temperature. We use the 
value of (a that results from the application of Gauss' 
principle of least constraint to the imposition of constant 
kinetic temperature: 



(a 



E£Li p? 



(22) 
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The atomic temperature Ta for a system of N a atoms 
with no internal degrees of freedom is defined as 



1 



(dN a - N c )k L 



E 

»=i 



(23) 



Here angled brackets denote an ensemble average, d is the 
dimensionality of the atomic system, N c is the number 
of constraints on the system (including constraints for 
conserved quantities). 

The equations of motion (EOM) for a molecular fluid 
can be formulated in either atomic or molecular ver- 
sions. In fact the molecular versions of the homogeneous 
isothermal EOM with a molecular thermostat at equilib- 
rium are similar to atomic EOM with a molecular ther- 
mostat, provided that all of the relevant forces are in- 
cluded [23]. The thermostatted EOM for molecular sys- 
tems are given by [29[ 



Pia 



Fm " Cm m" Pi 



(24) 



where ia denotes site a on molecule i. Ti a is the posi- 
tion, pi a is the momentum and m !Q is the mass of the 
site. The force on a site is separated into two terms: 
F ia is the contribution due to the Lennard- Jones type 
interactions on site a of molecule i and Ff a is either the 
constraint force or the bonding force. Cm is the molecular 
thermostat multiplier, given by 



Cm 



E^pIM 



(25) 



where Pj is the total force acting on molecule i and Mj is 
the mass of molecule i. Cm is derived from Gauss' prin- 
ciple of least constraint and acts to keep the molecular 
center of mass kinetic temperature Tm constant. Sev- 
eral algorithms are available for this purpose [47j]. The 
molecular temperature Tm is defined by 



T 



1 



M 



(dN m - N c )ki 



N 

E 

i=l 



(26) 



where N c is the number of translational center-of-mass 
degrees of freedom and depends on the total number 
of sites and the number of constraints on the system. 
The details of the constraint algorithm used to calcu- 
late Ff a have been discussed previously [H, [H, [46[ . All 
the systems were equilibrated for at least 10 6 time steps. 
The results from production runs were ensemble aver- 
aged over 14 runs, each of length 10 6 steps (i.e. a total 
of 1.4 x 10 7 time steps). The transverse momentum den- 
sity ACFs were computed over at least 20 reduced time 
units and the stress ACFs were computed over at least 
40 reduced time units. Both the transverse momentum 
density and stress ACFs were computed at wave-vectors 
k yn = lixnjLy where mode number n is from to 40 
with increment 2 and L y = [N a / p] 1 ^ 3 . For the remainder 



of this paper we drop the n index in k yn for simplic- 
ity. The ACFs were Laplace transformed with respect 
to time using Filon's rule |47|- We do not report the 
frequency dependent viscosities in this work. The wave- 
vector and frequency dependent viscosities were calcu- 
lated using Eqs. {T5| and (fTTjl . Eq. (|T3j) was used to 
obtain the non-zero wave-vector and frequency depen- 
dent viscosities and Eq. (| 1T[) was used to obtain the zero 
wave- vector viscosity. 

In this work, all quantities are expressed in reduced 
Lennard- Jones units. Thus, our reference units are 
the: reduced length r* = r/a, reduced number density 
p* = p/a 3 , reduced temperature T* = fc^T/e, reduced 
time t* = tj '(c(m/ 'e) 1 ^ 2 , reduced pressure P* = P(<7 3 /e), 
reduced energy E* = E/e and reduced viscosity r/* = 
T]o 2 1 Wjmej . For the remainder of this paper we apply 
these units and omit writing the asterisk. We will not 
distinguish between Tm and Ta, but simply use T to 
indicate the temperature. 



IV. RESULTS AND DISCUSSION 

The autocorrelation functions were evaluated for both 
108 and 864 molecule system in order to determine 
whether the results were system size dependent. There 
were no observed differences within the statistical errors 
in both monatomic and diatomic systems. We also note 
that in order to limit the number of figures, we do not 
display the results for the transverse momentum density 
and stress autocorrelation functions. However, we must 
mention that for monatomic systems both quantities (i.e. 
the transverse momentum density and stress ACFs) were 
in good agreement with those previously observed for 
Lennard- Jones monatomic liquids and their running in- 
tegrals have fully converged which suggests that the cor- 
relation functions have decayed to zero. 



A. Viscosity kernels in reciprocal space 

The reciprocal space kernels for monatomic and di- 
atomic fluids are plotted in figures ([T][3]). The error bars 
are smaller than the symbol sizes and therefore omit- 
ted in figures [2] and [3] Generally the statistical relia- 
bility of reciprocal space kernel data increases as k y in- 
creases. Our zero wave-vector, zero frequency viscosities 
for monatomic fluids agree well with those available in 
the literature. For the WCA system at the state point 
(p a = 0.375, T = 0.765) we found r? = 0.27±0.01 which 
agrees with the results of Hansen et al. 0.273 [12, l20l |. 
while at the state point (p a = 0.840, T = 1.0) we found 
r)o = 2.29±0.07, in agreement with the results of Matin et 
al. (2.1 ±0.2) For chlorine we found 770 = 6.89±0.32, 
which agrees with the limiting values (6.7±0.4) of the 
shear or elongational viscosities at zero strain rate ■ 

The wave-vector dependent viscosity for diatomic sys- 
tems, figure [TJ depicts a similar behaviour within the 
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TABLE II: Zero frequency, zero wave-vector shear viscosity and fitted parameter values for monatomic and diatomic systems 







WCA 


WCA WCA 


LJ LJ 


Chlorine 


State Point 


Pa 


0.375 


0.480 


U.o4U 


0.840 0.840 


1.088 




T 


0.765 


0.765 


1 nnn 


0.765 1.0 


0.97 


System size 


N a 




2048 






1728 




'/(> 


0.265(0.273[20j) 0.392 


9 9Qfl 

z.zyu 


2.810 2.650 


6.889 


2-term Gaussian, A2 = 1 - 


- Ai Eq. J33 A 


0.189(0.440[20]) 0.309 


U.100 


0.093 0.107 


0.407 






12.48(4.750[20J) 6.916 




10.04 9.088 


5.377 




£72 


2.1l6(1.376[20J) 


1.835 


9 C^09 


2.778 2.759 


1.236 




St- 


0.007(0.005[20J) 0.013 


U.U44 


0.021 0.027 


0.082 


2-term Gaussian Eq. (127[) 




0.792 


0.687 


U.O ( 4 


0.907 0.892 


0.592 




A, 


0.174 


0.254 


U.loo 


0.094 0.106 


0.407 




(7 1 


2.245 


2.113 


9 t;nn 

Z.Oyz 


2.776 2.765 


1.237 




il 2 


13.36 


7.745 


Q 1 9/1 


10.02 9.127 


5.377 




S r 


0.007 


0.011 


U.UoO 


0.022 0.031 


0.081 


4-term Gaussian Eq. (1271) 


A x 


0.432 


0.566 


0.778 


0.689 0.868 


0.398 




Ai 


0.394 


0.248 


0.118 


0.190 0.047 


0.538 




A 3 


0.120 


0.138 


0.088 


0.114 0.089 


0.055 




A 4 


0.056 


0.047 


0.017 


0.017 0.020 


0.008 




tr 1 


3.228 


2.826 


2.950 


2.709 2.814 


4.355 




<T 2 


1.261 


0.821 


0.651 


2.709 0.145 


1.155 




03 


8.165 


6.973 


8.496 


7.037 7.628 


10.46 




<T 4 


15.19 


14.74 


23.66 


24.94 19.99 


37.56 






0.008 


0.002 


0.012 


0.014 0.024 


0.018 


Lorentzian-type Eq. (128ft 


(1 


0.198(0.180[2pJ) 0.170 0.062 0.041 0.043 


0.239 






1.562(1.662[2pJ) 


1.715 


2.326 
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FIG. 1: f/(ky) versus k y for chlorine calculated using atomic 
and molecular formalisms (p a = 1.088, T = 0.97, N a = 1728). 



statistical uncertainty in both atomic and molecular for- 
malisms. 

It has been shown previously that numerous one pa- 
rameter functions failed to capture the behaviour of the 
reciprocal space kernel data [20] ■ We therefore present 
the best fits with two or more fitting parameters. We 
have identified two functional forms that fit the data well: 



an Nq term Gaussian function 

N G 

VG(k y ) = r)o A i c M-k 2 y /2a]) A h u j G M+ (27) 
3 

and a Lorentzian type function 

f ldk y ) = - "^—p a,/?eK + , (28) 

l + a|fcj,r 

We present the best fits of the data to (i) a two-term 
Gaussian function with freely estimated amplitudes (i.e. 
unconstrained fitting) termed as t)q 2 , (ii) to a two-term 
Gaussian function with interdependent amplitudes (i.e. 
constrained fitting \\^ G Aj = 1) given by Hansen et 
al. [2(| and termed as fjG 2 HJ ( m ) to a four-term Gaus- 
sian function with freely estimated amplitudes, termed as 
fjGi and (iv) to the Lorentzian type function, Eq. (|28|) . 
In order to measure the magnitude of the residuals we 
use the residual standard deviation defined as s r — 
\/Xm=i r2 / ( n s :~ n p ) where n s is the number of data 
points, n p is the number of fitting parameters, and r is 
the residual [48| . After an iterative curve fitting proce- 
dure the accurate estimation of rjo was kept fixed allowing 
all other parameters in Eqs. (|2T|) and ([28]) to be used as 
fitting parameters. In Table ITT1 we have listed the fitting 
parameters for monatomic and diatomic molecular fluids 
and compared to the previous results where possible. 
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FIG. 2: Normalized kernel data, best fit of Eq. |27)l with Ng = 2 and 4, and Eq. (|28|) and difference between the fits: (a) Best 
fits to normalized kernel for chlorine fluid (p a = 1.088, T = 0.97, N a = 1728); (b) Best fits to normalized kernel for LJ fluid 
(p a — 0.840, T = 1.0, N a — 2048); (c) Differences between the kernels and simulation data (a); (d) Differences between the 
fitted kernels (b). 



TABLE III: Total amplitude for Gaussian functional form 











WCA 




LJ 


Chlorine 


State Point 


pa 




0.375 


0.480 


0.840 


0.840 0.840 


1.088 




T 




0.765 


0.765 


1.000 


0.765 1.0 


0.97 


2-term Gaussian 


n= 




0.966 


0.941 


1.029 


1.001 0.998 


0.999 


4-term Gaussian 






1.002 


0.999 


1.001 


1.010 1.024 


0.999 



A useful check of the fitting can be performed by calcu- 
lating the total Gaussian amplitudes which should con- 
verge to the value of 1, Table fTTTl 

The reciprocal space results presented in figure for 
LJ and chlorine systems demonstrate that the four-term 
Gaussian function fits the data much better than the 
other two forms with a difference between the data and 
the fit of less than 0.5%, see figure [UJc). The two- 
term Gaussian fjG 2H fits the kernel data better than the 
Lorentzian-type function in the low-fc y region, figure[2ja), 
which suggests a more Gaussian-like behaviour in the 
low-fcy region, a fact previously observed by Hansen et. 
al [2fJ] for atomic fluids modeled with WCA potentials. 
Nevertheless the difference between the two-term Gaus- 
sian fit and data is less than 2% which still makes the 
fiG 2H a good analytical three parameter approximation 
of the reciprocal space viscosity kernel. The maximum 
difference between the Lorentzian-type fit and Gaussian 
fits are around 4% while the maximum difference between 
the Gaussian fits is about 2%, see figure [2Jd). Essen- 



tially, this suggests that, when computing the real space 
kernels, the four-term Gaussian functional form is to be 
trusted. It is obvious that eight parameters in the four- 
term Gaussian make its use less convenient, but on the 
other hand the Gaussian function can analytically be in- 
verse Fourier transformed while the inverse Fourier trans- 
form of the Lorentzian-type function can only be evalu- 
ated numerically for general values of /3. 

Figure EJa) shows the kernel data for a WCA fluid at 
two different densities along with two sets of data pub- 
lished previously by Hansen et al. [20l |. EMD is the 
set obtained from an equilibrium MD simulation at the 
same state point (p a = 0.375, T = 0.765) and NEMD is 
the set obtained from a nonequilibrium MD simulation 
based on the sinusoidal transverse force (STF) method. 
An excellent agreement between both sets of data was 
found. Figure [3]Jb) shows the normalized fit to Eq. (f2"5|). 
The normalized kernels, figure EJb), show a similar be- 
havior for k n < 4. Though the higher density kernel is 
slightly lower for k n > 4 they show a similar limiting be- 
havior. This effect was not seen by Hansen et al. due 
to lack of data for high wave vectors. Figure EJd) indi- 
cates that despite the difference between the interaction 
potentials, the results for LJ and WCA fluids are very 
close. This confirms that transport is dominated by re- 
pulsive interactions, rather than attractive. The sharper 
kernel for the diatomic system, figure EJd) , suggests a 
more Lorentzian-type behavior in the low wave-vector 
region. It is also important to mention that even though 
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FIG. 3: Reciprocal-space kernels of monatomic and diatomic fluids: (a) Kernel data of a WCA fluid at two different densities 
(T — 0.765, N a = 2048); (b) Best fit of normalized kernel data (a) to Lorentzian-type function Eq. (|28[): (c) Kernel data of 
a WCA fluid and LJ fluid at the same state point (p a = 0.840, T = 1.0, N a = 2048), Chlorine at (p a = 1.088, T = 0.97, 
N a = 1728); (d) Best fit of the normalized kernel data (c) to the Lorentzian-type function Eq. (|28[1 . 



the fitting parameters are significantly affected by tem- 
perature, the resulting kernels vary weakly over the range 
of temperatures chosen here. This was also observed by 
Hansen et al. for WCA monatomic fluids. 



B. Viscosity kernels in physical space 

The viscosity kernel in reciprocal space is an even func- 
tion since it is symmetric about the origin; thus the the 
real space kernel is symmetric because the Fourier trans- 
form keeps the even properties of the function. This 
means that the viscosity kernel in physical space can be 
found via an inverse Fourier cosine transform, F~ l [. . . ] (a 
special case of the continuous Fourier transform arising 
naturally when attempting to transform an even func- 
tion), of the viscosity kernel in reciprocal space. Since 
the integral is being computed over an interval symmet- 
ric about the origin (i.e. -oo to +oo), the second integral 
must vanish to zero, and the first may be simplified to 
give: 

oo 

F^iviky)} = V(y) = \/lJ fj(k y )cos(k y y)dk y . (29) 

o 

The inverse Fourier cosine transform of the Gaussian 
function, Eq. 1|2"7|). exists [4!| and it is even possible to 
obtain an analytical expression. For an Nq term Gaus- 



sian function the inverse Fourier cosine transform is 

N G 

m{v) = E A W cxp[-(a,2/) 2 /2] Aj, aj g R+. 

(30) 

Though the Lorentzian-type function given in Eq. (|28|) 
fulfills the criteria for having an inverse Fourier transform 
r]L(y) (i.e. the function is absolutely intcgrablc, square 
integrable and the function and its derivative are piece- 
wise continuous), the integral in Eq. (f2"9")) is not readily 
obtained analytically in the general case. However, the 
integral can be evaluated numerically. In this work, a 
Simpson method has been employed for this purpose. 

The real space kernels for atomic and diatomic flu- 
ids at zero frequency are presented in figure 0J [5] and |U 
Figure UJa) shows the resulting kernels for chlorine and 
figure IUb) shows the resulting kernels for a LJ fluid ex- 
tracted from two-term and four-term Gaussian functions 
Eq. ([30]) . and inverse Fourier transform of Eq. (f28|) . We 
find very little difference between the kernels obtained 
via two- and four-term Gaussians for these systems. Fig- 
ures EJc) and|l|d) show the differences between all three 
fits. It can be seen that there exists a significant differ- 
ence (almost 25%) between the kernels extracted from 
Gaussian and Lorentzian type functions for small y. The 
discrepancy decreases rapidly as y increases and becomes 
approximately zero for y > 1.5. The width of the kernel 
for chlorine is roughly 4-6 atomic diameters, figure EJa), 
and 3-5 atomic diameters for monatomic LJ and WCA 
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FIG. 4: Real space kernel of monatomic and diatomic fluids as predicted by Gaussian and Lorentzian fits of reciprocal space 
kernel data: (a) chlorine (p a = f .088, T — 0.97); (b) LJ (p a = 0.840, T = 1.0); (c) Differences between the kernels shown in 
(a); (d) Differences between the kernels shown in (b). 



fluids, figure @|b). 

For monatomic systems at relatively low densities 
(Pa = 0.375 — 0.480), the real space kernels are affected 
considerably by the functional form chosen to fit the re- 
ciprocal space kernel, figure [5) For instance, the equally 
weighted two-term Gaussian function, figure Oa), dis- 
torts the real space kernels and predicts a noticeably 
higher rj(y = 0) value. As we increase the density, the dis- 
crepancy between Gaussian functions, as well as between 
Gaussian and Lorcntzian-type functions, only partially 
reduces, figure EJb,d). The width of the kernel for WCA 
fluids at low density is roughly 2-4 atomic diameters. 

In figures Ela) and IH^b) we compare the unnormalized 
kernel data in y space extracted from four-term Gaussian 
and Lorentzian-type functional forms for all the simu- 
lated systems. Despite the fact that the difference be- 
tween the reciprocal kernels is less than 4%, figure [^d) 
(e.g. chlorine - dashed line), the kernels for the corre- 
sponding systems in real space look noticeably different 
(figurelUJa) and figurcEJb)), for all the systems (e.g. chlo- 
rine - dashed dotted line). However, the zero wave- vector 
viscosities obtained from both functional forms are very 
close, with less than 2% error. 

We can determine the zero wave- vector viscosities ?/o = 
r)(k = Q,uj = 0) by integrating the real space kernel over 
y, and thus test our numerical analysis. The zero wave- 
vector viscosity 770 obtained by a Gaussian function rjciy) 



TABLE IV: Effective viscosities evaluated from j?g 4 (j/), 
r/G 2H (y) and numerically from Tji(fe). The values can be 
compared to the zero frequency, zero wave-vector viscosities 
shown in Table ILT1 







WCA LJ 


Chlorine 


State Point 


Pa 


0.375 0.480 0.840 0.840 0.840 


1.088 




T 


0.765 0.765 1.000 0.765 1.000 


0.97 


2-term Gaussian 


>l« 


0.265 0.392 2.290 2.723 2.614 


6.881 


4-term Gaussian 


</<> 


0.265 0.390 2.288 2.807 2.653 


6.897 


Lorentzian 


1/0 


0.269 0.428 2.320 2.913 2.711 


7.049 



is 

00 00 
770= J r) G (y)dy=^jL J {Y^ A i a i ex P[-( a jy) 2 / 2 ]} d y- 

— OO —OO * 

(31) 

Since the general analytical expression for r]i,(y) does 
not exist [20j] we evaluate the integral numerically and 
present the results from all functional forms in Table llVl 
A comparison of the viscosities in Table ITVl with the sim- 
ulated zero frequency zero wave-vector shear viscosities 
given in Table HI] shows an integration error of less than 
3%. This confirms the accuracy of our numerical analysis 
techniques. 

It is of interest to discuss the real space viscosity ker- 
nels for monatomic and diatomic systems from a struc- 
tural point of view. For this purpose we define a struc- 
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FIG. 5: Real space kernel of monatomic WCA fluids as predicted by Gaussian fits of the reciprocal space kernel data, Eq. (|30|l : 
(a) Kernels obtained from two- and four-term Gaussian fits for a WCA fluid at p a = 0.375 and T = 0.765; (b) Kernels obtained 
from two- and four-term Gaussian fits for a WCA fluid at p a = 0.480 and T = 0.765; (c) Differences between the kernels for a 
WCA fluid at p a = 0.375 and T = 0.765; (d) Differences between the kernels for a WCA fluid at p a = 0.480 and T = 0.765. 
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FIG. 6: Real space viscosity kernels of monatomic and diatomic fluids, WCA (p a = 0.375, T = 0.765), WCA (p a — 0.480, 
T = 0.765), WCA (p a = 0.840, T = 1.0), LJ (p a = 0.840, T = 1.0), chlorine (p a = 1.088, T = 0.97): (a) Kernels obtained from 
the four-term Gaussian functional form Eq. <)30[) ; (b) Kernels obtained numerically from the Lorentzian-type functional form 
Eq. ([281). 



tural normalization factor 

oo 

/ r l9( r ) ~ l] 2 dr 

a a = ^ (32) 

J[ 9 (r)-l?dr 
o 

where g(r) is the radial distribution function (RDF). 
Eq. (|32|) is a measure of the range over which the cor- 
relation function decays to zero and therefore could be 
regarded as a correlation length of the radial distribu- 
tion function. The RDF (or structure factor in reciprocal 



space) can be defined either in terms of the vector norm 
rij between the atoms i and j or between the centres of 

mass of molecules i and j: g(r) = ^ ^'= 1 ^j>i^ r r 'jD ^ 

where N is the total number of atoms or molecules, and 
p is the atomic or molecular number density. 

The radial distribution functions and normalization 
factors are presented in figure [7| We can see that the 
RDF are typical monatomic and diatomic Lennard- Jones 
pair correlation functions. £ ff generally increases as we in- 
crease the density and temperature from 0.605 at state 
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FIG. 7: Pair-distance correlation function g(r) and normaliza- 
tion factors £ g , Eq. ([32]): WCA [(a) p a = 0.375, (b) p a = 0.480 
both at T = 0.765], (c) WCA (p a = 0.840, T = 1.0); LJ (d) 
(p a = 0.840, T = 1.0); chlorine (e) (p a = 1.088, T = 0.97). 
For clarity, the RDFs are shifted upwards by 3 units. 



point p a = 0.375, T = 0.756 to 0.730 at p a = 0.480, 
T = 1.0 and only slightly increases as we increase the 
cutoff distance, i.e. switch from WCA system to a LJ 
system at the same state point. £ g for chlorine at state 
point p a = 1.088, T = 0.97 was found to be 0.585. 

The normalized kernels with respect to rj(y = 0) and 
normalization factor £ g are shown in figures [5Ja) and 
[E^b). While generally the width of unnormalized ker- 
nels increases as we increase the density (figures Eta) and 
|6jb)), the width of the normalized kernels of WCA flu- 
ids decreases marginally as we increase the density from 
0.376 (continuous line) to 0.840 (short-dashed line). The 
LJ system shows a slightly narrower kernel (dotted line) 
compared to the WCA system at the same state point 
(short-dashed line). Though the kernels obtained from 
both functional forms are quite close to each other (al- 
most identical for values of y of about half of the atomic 
or molecular diameters, i.e. y = 0.5er), wc can see in 
figure Hlja) and [5Jb) that the structural normalization 
did not completely remove the discrepancy between the 
normalized kernels of the WCA system at different densi- 
ties, and the normalized kernels of WCA, LJ and chlorine 
systems, for values higher than y = a . If wc recall that 
figure [E^b) is based on a Lorcntzian-type fit, a further 
question as to whether the kernel differences are due to 
numerical analysis, i.e. the choice of the fitting function 



or due to improper structural factor arises. A four-term 
Gaussian only shows a slightly narrower kernels which 
suggests a need for a more comprehensive structural nor- 
malization. 



V. CONCLUSION 

The wave-vector dependent viscosity of monatomic 
Lennard- Jones, monatomic Weeks-Chandler-Andersen 
and diatomic (liquid chlorine) fluids over a large wave- 
vector range and for a variety of state points has been 
computed. The equilibrium molecular dynamics calcula- 
tion involved the evaluation of the transverse momentum 
density and shear stress autocorrelation functions in both 
atomic and molecular hydrodynamic representations for 
molecular fluids. The main results can be summarized as 
follows: 

(i) For monatomic fluids the shape of the normalized 
viscosity kernel in reciprocal space in the low wave- vector 
region is the same for a whole range of densities consid- 
ered here. Though the normalized reciprocal kernels in- 
significantly decreases with the density they show a sim- 
ilar limiting behaviour at high k y values. For the LJ 
potential compared to a WCA potential we find higher 
viscosities in the low wave-vector region but the normal- 
ized shape of the kernels are almost identical. 

(ii) For liquid chlorine, the wave-vector dependent vis- 
cosity shows a similar behaviour in both atomic and 
molecular formalisms within statistical uncertainty. 

(iii) While a relatively simple Lorcntzian-type function 
fits the atomic and diatomic data well over the entire 
range of k y at all the state points it is not possible to 
analytically inverse Fourier transform it to the real space 
domain. Therefore one may consider an expansion up 
to a four-term Gaussian which gives better accuracy in 
reciprocal space compared to the Lorentzian-type func- 
tion. Our analysis of the high k y regime reveals that the 
two-term equally weighted Gaussian functional form is 
inaccurate in predicting the real space kernels whilst the 
unequally weighted Gaussian only slightly improves the 
fit. 

(iv) The overall conclusion is that the real space viscos- 
ity kernel for chlorine has a width of roughly 4-6 atomic 
diameters while for monatomic systems at high densities 
the width is about 3-5 atomic diameters and 2-4 atomic 
diameters at low densities. This means that generalized 
hydrodynamics must be used in predicting the flow prop- 
erties of molecular fluids on length scales where the gra- 
dient in the strain rate varies significantly on these scales. 
Consequently a nonlocal constitutive equations should be 
invoked for a complete description of flows at atomic and 
molecular scales under such conditions. 

Finally, our results for molecular fluids should also pro- 
vide a good test for more complex molecular systems and 
the methodology can easily be used for instance in chain- 
like molecules. 
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FIG. 8: Normalized real space viscosity kernels of monatomic and diatomic fluids, WCA (p a = 0.375, T = 0.765), WCA 
(p a = 0.480, T = 0.765), WCA (p a = 0.840, T = 1.0), LJ (p a = 0.840, T = 1.0), chlorine (p a = 1.088, T = 0.97): 
(a) Normalized kernels, shown in Fig. HJa), obtained from the four-term Gaussian functional form Eq. (|30p : (b) Normalized 
kernels, shown in Fig. |6jb), obtained numerically from the Lorentzian-type functional form Eq. I)28[l. 
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